Introduction
Packages
library(tidyverse)
library(plotly)
library(scales)
library(gganimate)
library(knitr)
library(DT)
library(colorRamps)Input Parameters
S <- 1370
A <- 204
B <- 2.17
K <- 3.86
ai <- 0.62
ab <- 0.25
aW <- 0.75
aB <- 0.25
gamma <- 2.2
delta <- 10/gamma
c <- 7
k <- 0.003265*0.75
T0 <- 20
D <- 0.3 # Death Rate
w0 <- 0.5
b0 <- 0.2
u0 <- 1-w0-b0Fit from Textbook Spreadsheet
zones <- c(85,75,65,55,45,35,25,15,6)
coszones <- cospi(zones/180)
sunWt <- c(0.5,0.531,0.624,0.77,0.892,1.021,1.12,1.189,1.219)
df <- data.frame(zones,sunWt)
f <- function(zones,d,n,k){
d*cos(n*zones)^2+k } #cos or cos^2 ??
Fit <- nls(sunWt ~ f(zones,d,n,k),data= df, start=list(d=1,n=0.015,k=0.3))
ggplot(df, aes(zones,sunWt))+geom_point()+geom_smooth(aes(y=f(zones,0.7768699,0.0164348,0.4617747)))+xlab("Latitude")+ylab("Solar Weight Factor")
Functional Forms
gauss <- function(x,m,sd,b){
((24+b)/(0.00798*sqrt(2*pi*sd^2)))*exp(-(x-m)^2/(2*sd^2))-b
}
#SOLAR LUMINOSITY (Latitude)
Func <- function(x){ 0.7768699*cos(0.0164348*x)^2+0.4617747 }
# SOLAR LUMINOSITY (Time)
Sun1 <- function(x,a){a*x/100}
Sun2 <- function(x){1370*(sinpi((x+90)/180))^2}
Sun3 <- function(x){1370-(((1370)/(sqrt(2*pi*1^2)))*exp(-(x-50)^2/(2*1^2)))}
Sun4 <- function(x){ifelse(x==50,1370/3,1370)}
Sun5 <- function(x){(1/100) * (abs(x-150)+25)}
Sun6 <- function(x){1370*(1+0.1*cospi(x/180))}
Incident <- function(x,y){x*y/4}
# ALBEDO
Step <- function(x,c){ifelse(x<c, 0.6, 0.3)}
alb <- function(x,a,b,c,d){ (exp(c*(x+d)) / (exp(c*(x+d))+1)) * (b-a) + a }
# EBM
ebm01 <- function(cycles,A,B,K,ai,ab,gamma,delta) {
Incident <- function(x,y){ x*y/4 }
Func <- function(x){ 0.7768699*cos(0.0164348*x)^2+0.4617747 }
Zones <- seq(-89, 89, by = 2)
cosZones <- abs(cospi(Zones/180))
t <- c(1:cycles)
Temperature <- rep(0,cycles)
Ti <- gauss(Zones,0,50,31.6)
SunWt <- Func(Zones)
Rin <- Incident(S,SunWt)
T <- Ti
a <- alb(T,ai,ab,gamma,delta)
for(i in t) { Tcos <- cosZones*T
Tm <- sum(Tcos)/sum(cosZones)
T <- (Rin*(1-a)+K*Tm-A) / (B+K)
a <- alb(T,ai,ab,gamma,delta)
Temperature[i] <- Tm }
return( data.frame(Zones,T,a,Ti)) }
ebm02 <- function(cycles,A,B,K,ai,ab,gamma,delta) {
Incident <- function(x,y){ x*y/4 }
Func <- function(x){ 0.7768699*cos(0.0164348*x)^2+0.4617747 }
Zones <- seq(-89, 89, by = 2)
cosZones <- abs(cospi(Zones/180))
t <- c(1:cycles)
Temperature <- rep(0,cycles)
Ti <- gauss(Zones,0,50,31.6)
SunWt <- Func(Zones)
Rin <- Incident(S,SunWt)
T <- Ti
a <- alb(T,ai,ab,gamma,delta)
for(i in t)
{Tcos <- cosZones*T
Tm <- sum(Tcos)/sum(cosZones)
T <- (Rin*(1-a)+K*Tm-A) / (B+K)
a <- alb(T,ai,ab,gamma,delta)
Temperature[i] <- Tm }
return( data.frame(t,Temperature) ) }
ebm11 <- function(cycles1,cycles2,A,B,K,ai,ab,gamma,delta) {
Incident <- function(x,y){ x*y/4 }
Func <- function(x){ 0.7768699*cos(0.0164348*x)^2+0.4617747 }
Sun1 <- function(x,a){a*x/100}
Zones <- seq(-89, 89, by = 2)
cosZones <- abs(cospi(Zones/180))
J <- rep(0,cycles1)
TEMP1 <- matrix(NA, nrow=90, ncol=cycles1)
Temp <- rep(0,cycles1)
Sun1 <- function(x,a){a*x/100}
SunWt <- Func(Zones)
for(j in c(1:cycles1)){
T <- gauss(Zones,0,50,31.6)
a <- alb(T,ai,ab,gamma,delta)
S <- Sun1(1370,j)
Rin <- Incident(S,SunWt)
for(i in c(1:cycles2))
{Tcos <- cosZones*T
Tm <- sum(Tcos)/sum(cosZones)
T <- (Rin*(1-a)+K*Tm-A) / (B+K)
a <- alb(T,ai,ab,gamma,delta)
TEMP1[,j] <- T
J[j] <- j
Temp[j] <- Tm } }
return( data.frame(J,Temp) ) }
ebm12 <- function(cycles1,cycles2,A,B,K,ai,ab,gamma,delta) {
Incident <- function(x,y){ x*y/4 }
Func <- function(x){ 0.7768699*cos(0.0164348*x)^2+0.4617747 }
Sun1 <- function(x,a){a*x/100}
Zones <- seq(-89, 89, by = 2)
cosZones <- abs(cospi(Zones/180))
SunWt <- Func(Zones)
J <- rep(0,cycles1)
TEMP1 <- matrix(NA, nrow=90, ncol=cycles1)
Temp <- rep(0,cycles1)
for(j in c(1:cycles1)){
T <- gauss(Zones,0,50,31.6)
a <- alb(T,ai,ab,gamma,delta)
S <- Sun1(1370,j)
Rin <- Incident(S,SunWt)
for(i in c(1:cycles2))
{Tcos <- cosZones*T
Tm <- sum(Tcos)/sum(cosZones)
T <- (Rin*(1-a)+K*Tm-A) / (B+K)
a <- alb(T,ai,ab,gamma,delta)
TEMP1[,j] <- T
J[j] <- j
Temp[j] <- Tm } }
return( TEMP1 ) }
ebm21 <- function(cycles1,cycles2,A,B,K,ai,ab,gamma,delta ) {
Incident <- function(x,y){ x*y/4 }
Func <- function(x){ 0.7768699*cos(0.0164348*x)^2+0.4617747 }
Sun2 <- function(x){1370*(sinpi((x+90)/180))^2}
Zones <- seq(-89, 89, by = 2)
cosZones <- abs(cospi(Zones/180))
SunWt <- Func(Zones)
J <- rep(0,cycles1)
Temp2 <- rep(0,cycles1)
T2 <- gauss(Zones,0,50,31.6)
TEMP2 <- matrix(NA, nrow=90, ncol=cycles1)
a <- alb(T2,ai,ab,gamma,delta)
Sarr <- rep(0,cycles1)
for(j in c(1:cycles1)){
S <- Sun2(j)
Rin <- Incident(S,SunWt)
for(i in c(1:cycles2))
{Tcos <- cosZones*T2
Tm <- sum(Tcos)/sum(cosZones)
T2 <- (Rin*(1-a)+K*Tm-A) / (B+K)
a <- alb(T2,ai,ab,gamma,delta)
}
Sarr[j] <- Sun2(j)
TEMP2[,j] <- T2
J[j] <- j
Temp2[j] <- Tm }
return( data.frame(J,Temp2) ) }
ebm22 <- function(cycles1,cycles2,A,B,K,ai,ab,gamma,delta ) {
Incident <- function(x,y){ x*y/4 }
Func <- function(x){ 0.7768699*cos(0.0164348*x)^2+0.4617747 }
Sun2 <- function(x){1370*(sinpi((x+90)/180))^2}
Zones <- seq(-89, 89, by = 2)
cosZones <- abs(cospi(Zones/180))
J <- rep(0,cycles1)
Temp2 <- rep(0,cycles1)
T2 <- gauss(Zones,0,50,31.6)
TEMP2 <- matrix(NA, nrow=90, ncol=cycles1)
a <- alb(T2,ai,ab,gamma,delta)
Sarr <- rep(0,cycles1)
SunWt <- Func(Zones)
for(j in c(1:cycles1)){
S <- Sun2(j)
Rin <- Incident(S,SunWt)
for(i in c(1:cycles2))
{Tcos <- cosZones*T2
Tm <- sum(Tcos)/sum(cosZones)
T2 <- (Rin*(1-a)+K*Tm-A) / (B+K)
a <- alb(T2,ai,ab,gamma,delta)
}
Sarr[j] <- Sun2(j)
TEMP2[,j] <- T2
J[j] <- j
Temp2[j] <- Tm }
return( TEMP2 ) }
ebm23 <- function(cycles1,cycles2,A,B,K,ai,ab,gamma,delta ) {
Incident <- function(x,y){ x*y/4 }
Func <- function(x){ 0.7768699*cos(0.0164348*x)^2+0.4617747 }
Sun2 <- function(x){1370*(sinpi((x+90)/180))^2}
Zones <- seq(-89, 89, by = 2)
cosZones <- abs(cospi(Zones/180))
J <- rep(0,cycles1)
Temp2 <- rep(0,cycles1)
T2 <- gauss(Zones,0,50,31.6)
TEMP2 <- matrix(NA, nrow=90, ncol=cycles1)
a <- alb(T2,ai,ab,gamma,delta)
Sarr <- rep(0,cycles1)
SunWt <- Func(Zones)
for(j in c(1:cycles1)){
S <- Sun2(j)
Rin <- Incident(S,SunWt)
for(i in c(1:cycles2))
{Tcos <- cosZones*T2
Tm <- sum(Tcos)/sum(cosZones)
T2 <- (Rin*(1-a)+K*Tm-A) / (B+K)
a <- alb(T2,ai,ab,gamma,delta)
}
Sarr[j] <- Sun2(j)
TEMP2[,j] <- T2
J[j] <- j
Temp2[j] <- Tm }
return( data.frame(Sarr,Temp2) ) }
ebm31 <- function(cycles1,cycles2,A,B,K,ai,ab,gamma,delta ) {
Incident <- function(x,y){ x*y/4 }
Func <- function(x){ 0.7768699*cos(0.0164348*x)^2+0.4617747 }
Sun3 <- function(x){1370-(((1370)/(sqrt(2*pi*1^2)))*exp(-(x-50)^2/(2*1^2)))}
Zones <- seq(-89, 89, by = 2)
cosZones <- abs(cospi(Zones/180))
SunWt <- Func(Zones)
J <- rep(0,cycles1)
TEMP3 <- matrix(NA, nrow=90, ncol=cycles1)
Temp <- rep(0,cycles1)
T3 <- gauss(Zones,0,50,31.6)
a <- alb(T3,ai,ab,gamma,delta)
for(j in c(1:cycles1)){
S <- Sun3(j)
Rin <- Incident(S,SunWt)
for(i in c(1:cycles2))
{Tcos <- cosZones*T3
Tm <- sum(Tcos)/sum(cosZones)
T3 <- (Rin*(1-a)+K*Tm-A) / (B+K)
a <- alb(T3,ai,ab,gamma,delta)
}
TEMP3[,j] <- T3
J[j] <- j
Temp[j] <- Tm
}
return( data.frame(J,Temp) ) }
ebm32 <- function(cycles1,cycles2,A,B,K,ai,ab,gamma,delta ) {
Incident <- function(x,y){ x*y/4 }
Func <- function(x){ 0.7768699*cos(0.0164348*x)^2+0.4617747 }
Sun3 <- function(x){1370-(((1370)/(sqrt(2*pi*1^2)))*exp(-(x-50)^2/(2*1^2)))}
Zones <- seq(-89, 89, by = 2)
cosZones <- abs(cospi(Zones/180))
SunWt <- Func(Zones)
J <- rep(0,cycles1)
TEMP3 <- matrix(NA, nrow=90, ncol=cycles1)
Temp <- rep(0,cycles1)
T3 <- gauss(Zones,0,50,31.6)
a <- alb(T3,ai,ab,gamma,delta)
for(j in c(1:cycles1)){
S <- Sun3(j)
Rin <- Incident(S,SunWt)
for(i in c(1:cycles2))
{Tcos <- cosZones*T3
Tm <- sum(Tcos)/sum(cosZones)
T3 <- (Rin*(1-a)+K*Tm-A) / (B+K)
a <- alb(T3,ai,ab,gamma,delta)
}
TEMP3[,j] <- T3
J[j] <- j
Temp[j] <- Tm
}
return( TEMP3 ) }
ebm41 <- function(cycles1,cycles2,A,B,K,ai,ab,gamma,delta ) {
Incident <- function(x,y){ x*y/4 }
Func <- function(x){ 0.7768699*cos(0.0164348*x)^2+0.4617747 }
Sun4 <- function(x){ifelse(x==50,1370/3,1370)}
Zones <- seq(-89, 89, by = 2)
cosZones <- abs(cospi(Zones/180))
SunWt <- Func(Zones)
J <- rep(0,cycles1)
Temp <- rep(0,cycles1)
T4 <- gauss(Zones,0,50,31.6)
TEMP4 <- matrix(NA, nrow=90, ncol=cycles1)
a <- alb(T4,ai,ab,gamma,delta)
for(j in c(1:cycles1)){
S <- Sun4(j)
Rin <- Incident(S,SunWt)
for(i in c(1:cycles2))
{Tcos <- cosZones*T4
Tm <- sum(Tcos)/sum(cosZones)
T4 <- (Rin*(1-a)+K*Tm-A) / (B+K)
a <- alb(T4,ai,ab,gamma,delta)
}
TEMP4[,j] <- T4
J[j] <- j
Temp[j] <- Tm
}
return( data.frame(J,Temp) ) }
ebm42 <- function(cycles1,cycles2,A,B,K,ai,ab,gamma,delta ) {
Incident <- function(x,y){ x*y/4 }
Func <- function(x){ 0.7768699*cos(0.0164348*x)^2+0.4617747 }
Sun4 <- function(x){ifelse(x==50,1370/3,1370)}
Zones <- seq(-89, 89, by = 2)
cosZones <- abs(cospi(Zones/180))
SunWt <- Func(Zones)
J <- rep(0,cycles1)
Temp <- rep(0,cycles1)
T4 <- gauss(Zones,0,50,31.6)
TEMP4 <- matrix(NA, nrow=90, ncol=cycles1)
a <- alb(T4,ai,ab,gamma,delta)
for(j in c(1:cycles1)){
S <- Sun4(j)
Rin <- Incident(S,SunWt)
for(i in c(1:cycles2))
{Tcos <- cosZones*T4
Tm <- sum(Tcos)/sum(cosZones)
T4 <- (Rin*(1-a)+K*Tm-A) / (B+K)
a <- alb(T4,ai,ab,gamma,delta)
}
TEMP4[,j] <- T4
J[j] <- j
Temp[j] <- Tm
}
return( TEMP4 ) }
ebm51 <- function(cycles1,cycles2,A,B,K,ai,ab,gamma,delta ) {
Incident <- function(x,y){ x*y/4 }
Func <- function(x){ 0.7768699*cos(0.0164348*x)^2+0.4617747 }
Sun5 <- function(x){(1/100) * (abs(x-150)+25)}
Zones <- seq(-89, 89, by = 2)
cosZones <- abs(cospi(Zones/180))
SunWt <- Func(Zones)
J <- rep(0,cycles1)
Temp5 <- rep(0,cycles1)
T5 <- gauss(Zones,0,50,31.6)
TEMP5 <- matrix(NA, nrow=90, ncol=cycles1)
a <- alb(T5,ai,ab,gamma,delta)
Sarr <- rep(0,cycles1)
for(j in c(0:cycles1)){
S <- Sun5(j)
Rin <- 1370*Incident(S,SunWt)
for(i in c(1:cycles2))
{Tcos <- cosZones*T5
Tm <- sum(Tcos)/sum(cosZones)
T5 <- (Rin*(1-a)+K*Tm-A) / (B+K)
a <- alb(T5,ai,ab,gamma,delta)
}
Sarr[j] <- Sun5(j)
J[j] <- j
Temp5[j] <- Tm
}
return( data.frame(Sarr,Temp5) ) }
ebm_ND1 <- function(cycles,w0,b0,A,B,K,ai,ab,aW,aB,gamma,delta){
Incident <- function(x,y){ x*y/4 }
Func <- function(x){ 0.7768699*cos(0.0164348*x)^2+0.4617747 }
Sun6 <- function(x){1370*(1+0.1*cospi(x/180))}
Zones <- seq(-89, 89, by = 2)
cosZones <- abs(cospi(Zones/180))
SunWt <- Func(Zones)
Rin <- Incident(S,SunWt)
Sun6 <- function(x){1370*(1+0.1*cospi(x/180))}
T <- gauss(Zones,0,50,31.6)-6
w <- rep(w0,length(Zones)) #0.5
b <- rep(b0,length(Zones)) #0.2
u <- rep(1-w0-b0,length(Zones))
a <- w*aW+b*aB+u*alb(T,ai,ab,gamma,delta)
Barr <- rep(0,cycles)
Warr <- rep(0,cycles)
Uarr <- rep(0,cycles)
Tarr <- rep(0,cycles)
I <- rep(0,cycles)
TEMP <- matrix(NA, nrow=90, ncol=cycles)
for(i in c(1:cycles)) {
S <- Sun6(i) # oppure costante S <- 1370
Rin <- Incident(S,SunWt)
Tcos <- cosZones*T
Tm <- sum(Tcos)/sum(cosZones)
T <- (Rin*(1-a)+K*Tm-A) / (B+K)
TEMP[,i] <- T
a <- alb(T,ai,ab,gamma,delta)
I[i] <- i
Tarr[i] <- T[45]
}
return( data.frame(Zones,w,b,u,T) )}
ebm_ND2 <- function(cycles,w0,b0,A,B,K,ai,ab,aW,aB,gamma,delta){
Incident <- function(x,y){ x*y/4 }
Func <- function(x){ 0.7768699*cos(0.0164348*x)^2+0.4617747 }
Sun6 <- function(x){1370*(1+0.1*cospi(x/180))}
Zones <- seq(-89, 89, by = 2)
cosZones <- abs(cospi(Zones/180))
SunWt <- Func(Zones)
Rin <- Incident(S,SunWt)
Sun6 <- function(x){1370*(1+0.1*cospi(x/180))}
T <- gauss(Zones,0,50,31.6)-6
w <- rep(w0,length(Zones)) #0.5
b <- rep(b0,length(Zones)) #0.2
u <- rep(1-w0-b0,length(Zones))
a <- w*aW+b*aB+u*alb(T,ai,ab,gamma,delta)
Barr <- rep(0,cycles)
Warr <- rep(0,cycles)
Uarr <- rep(0,cycles)
Tarr <- rep(0,cycles)
I <- rep(0,cycles)
TEMP <- matrix(NA, nrow=90, ncol=cycles)
for(i in c(1:cycles)) {
S <- Sun6(i) # oppure costante S <- 1370
Rin <- Incident(S,SunWt)
Tcos <- cosZones*T
Tm <- sum(Tcos)/sum(cosZones)
T <- (Rin*(1-a)+K*Tm-A) / (B+K)
TEMP[,i] <- T
a <- alb(T,ai,ab,gamma,delta)
I[i] <- i
Tarr[i] <- T[45]
}
return( data.frame(I,Barr,Warr,Uarr,Tarr) )}
ebm_ND3 <- function(cycles,w0,b0,A,B,K,ai,ab,aW,aB,gamma,delta){
Incident <- function(x,y){ x*y/4 }
Func <- function(x){ 0.7768699*cos(0.0164348*x)^2+0.4617747 }
Sun6 <- function(x){1370*(1+0.1*cospi(x/180))}
Zones <- seq(-89, 89, by = 2)
cosZones <- abs(cospi(Zones/180))
SunWt <- Func(Zones)
Rin <- Incident(S,SunWt)
Sun6 <- function(x){1370*(1+0.1*cospi(x/180))}
T <- gauss(Zones,0,50,31.6)-6
w <- rep(w0,length(Zones)) #0.5
b <- rep(b0,length(Zones)) #0.2
u <- rep(1-w0-b0,length(Zones))
a <- w*aW+b*aB+u*alb(T,ai,ab,gamma,delta)
Barr <- rep(0,cycles)
Warr <- rep(0,cycles)
Uarr <- rep(0,cycles)
Tarr <- rep(0,cycles)
I <- rep(0,cycles)
TEMP <- matrix(NA, nrow=90, ncol=cycles)
for(i in c(1:cycles)) {
S <- Sun6(i) # oppure costante S <- 1370
Rin <- Incident(S,SunWt)
Tcos <- cosZones*T
Tm <- sum(Tcos)/sum(cosZones)
T <- (Rin*(1-a)+K*Tm-A) / (B+K)
TEMP[,i] <- T
a <- alb(T,ai,ab,gamma,delta)
I[i] <- i
Tarr[i] <- T[45]
}
return( TEMP )}
ebm_D1 <- function(cycles,w0,b0,c,k,D,A,B,K,ai,ab,aW,aB,gamma,delta){
Incident <- function(x,y){ x*y/4 }
Func <- function(x){ 0.7768699*cos(0.0164348*x)^2+0.4617747 }
Sun6 <- function(x){1370*(1+0.1*cospi(x/180))}
Zones <- seq(-89, 89, by = 2)
cosZones <- abs(cospi(Zones/180))
SunWt <- Func(Zones)
Rin <- Incident(S,SunWt)
T <- gauss(Zones,0,50,31.6)-6
w <- rep(w0,length(Zones)) #0.5
b <- rep(b0,length(Zones)) #0.2
u <- rep(1-w0-b0,length(Zones))
a <- w*aW+b*aB+u*alb(T,ai,ab,gamma,delta)
Barr <- rep(0,cycles)
Warr <- rep(0,cycles)
Uarr <- rep(0,cycles)
Tarr <- rep(0,cycles)
I <- rep(0,cycles)
TEMP <- matrix(NA, nrow=length(Zones), ncol=cycles)
for(i in c(1:cycles)) {
S <- Sun6(i) # oppure costante S <- 1370
Rin <- Incident(S,SunWt)
Tcos <- cosZones*T
Tm <- sum(Tcos)/sum(cosZones)
T <- (Rin*(1-a)+K*Tm-A) / (B+K)
TEMP[,i] <- T
Tw <- T+c*(a-aW)
Tb <- T+c*(a-aB)
Fw <- 1-k*(T0-Tw)^2
Fb <- 1-k*(T0-Tb)^2
for(j in c(1:length(Zones))){
if(Fw[j]<0){Fw[j]=0}
if(Fb[j]<0){Fb[j]=0} }
w <- w+w*(u*Fw-D)
b <- b+b*(u*Fb-D)
for(j in c(1:length(Zones))){
if(w[j]<0.001){w[j]=0.001}
if(b[j]<0.001){b[j]=0.001} }
u <- 1-w-b
a <- w*aW+b*aB+u*alb(T,ai,ab,gamma,delta)
Barr[i] <- b[45]
Warr[i] <- w[45]
Uarr[i] <- u[45]
I[i] <- i
Tarr[i] <- T[45]
}
return ( data.frame(Zones,w,b,u,T) )
}
ebm_D2 <- function(cycles,w0,b0,c,k,D,A,B,K,ai,ab,aW,aB,gamma,delta){
Incident <- function(x,y){ x*y/4 }
Func <- function(x){ 0.7768699*cos(0.0164348*x)^2+0.4617747 }
Sun6 <- function(x){1370*(1+0.1*cospi(x/180))}
Zones <- seq(-89, 89, by = 2)
cosZones <- abs(cospi(Zones/180))
SunWt <- Func(Zones)
Rin <- Incident(S,SunWt)
Sun6 <- function(x){1370*(1+0.1*cospi(x/180))}
T <- gauss(Zones,0,50,31.6)-6
w <- rep(w0,length(Zones)) #0.5
b <- rep(b0,length(Zones)) #0.2
u <- rep(1-w0-b0,length(Zones))
a <- w*aW+b*aB+u*alb(T,ai,ab,gamma,delta)
Barr <- rep(0,cycles)
Warr <- rep(0,cycles)
Uarr <- rep(0,cycles)
Tarr <- rep(0,cycles)
I <- rep(0,cycles)
TEMP <- matrix(NA, nrow=90, ncol=cycles)
for(i in c(1:cycles)) {
S <- Sun6(i) # oppure costante S <- 1370
Rin <- Incident(S,SunWt)
Tcos <- cosZones*T
Tm <- sum(Tcos)/sum(cosZones)
T <- (Rin*(1-a)+K*Tm-A) / (B+K)
TEMP[,i] <- T
Tw <- T+c*(a-aW)
Tb <- T+c*(a-aB)
Fw <- 1-k*(T0-Tw)^2
Fb <- 1-k*(T0-Tb)^2
for(j in c(1:length(Zones))){
if(Fw[j]<0){Fw[j]=0}
if(Fb[j]<0){Fb[j]=0} }
w <- w+w*(u*Fw-D)
b <- b+b*(u*Fb-D)
for(j in c(1:length(Zones))){
if(w[j]<0.001){w[j]=0.001}
if(b[j]<0.001){b[j]=0.001} }
u <- 1-w-b
a <- w*aW+b*aB+u*alb(T,ai,ab,gamma,delta)
Barr[i] <- b[45]
Warr[i] <- w[45]
Uarr[i] <- u[45]
I[i] <- i
Tarr[i] <- T[45]
}
return ( data.frame(I,Barr,Warr,Uarr,Tarr) )
}
ebm_D3 <- function(cycles,w0,b0,c,k,D,A,B,K,ai,ab,aW,aB,gamma,delta){
Incident <- function(x,y){ x*y/4 }
Func <- function(x){ 0.7768699*cos(0.0164348*x)^2+0.4617747 }
Sun6 <- function(x){1370*(1+0.1*cospi(x/180))}
Zones <- seq(-89, 89, by = 2)
cosZones <- abs(cospi(Zones/180))
SunWt <- Func(Zones)
Rin <- Incident(S,SunWt)
Sun6 <- function(x){1370*(1+0.1*cospi(x/180))}
T <- gauss(Zones,0,50,31.6)-6
w <- rep(w0,length(Zones)) #0.5
b <- rep(b0,length(Zones)) #0.2
u <- rep(1-w0-b0,length(Zones))
a <- w*aW+b*aB+u*alb(T,ai,ab,gamma,delta)
Barr <- rep(0,cycles)
Warr <- rep(0,cycles)
Uarr <- rep(0,cycles)
Tarr <- rep(0,cycles)
I <- rep(0,cycles)
TEMP <- matrix(NA, nrow=length(Zones), ncol=cycles)
for(i in c(1:cycles)) {
S <- Sun6(i) # oppure costante S <- 1370
Rin <- Incident(S,SunWt)
Tcos <- cosZones*T
Tm <- sum(Tcos)/sum(cosZones)
T <- (Rin*(1-a)+K*Tm-A) / (B+K)
TEMP[,i] <- T
Tw <- T+c*(a-aW)
Tb <- T+c*(a-aB)
Fw <- 1-k*(T0-Tw)^2
Fb <- 1-k*(T0-Tb)^2
for(j in c(1:length(Zones))){
if(Fw[j]<0){Fw[j]=0}
if(Fb[j]<0){Fb[j]=0} }
w <- w+w*(u*Fw-D)
b <- b+b*(u*Fb-D)
for(j in c(1:length(Zones))){
if(w[j]<0.001){w[j]=0.001}
if(b[j]<0.001){b[j]=0.001} }
u <- 1-w-b
a <- w*aW+b*aB+u*alb(T,ai,ab,gamma,delta)
Barr[i] <- b[45]
Warr[i] <- w[45]
Uarr[i] <- u[45]
I[i] <- i
Tarr[i] <- T[45]
}
return ( TEMP )
}
ebm_Db1 <- function(cycles,s,w0,b0,c,k,D,A,B,K,ai,ab,aW,aB,gamma,delta){
Incident <- function(x,y){ x*y/4 }
Func <- function(x){ 0.7768699*cos(0.0164348*x)^2+0.4617747 }
Sun7 <- function(x,y){x*y}
Zones <- seq(-89, 89, by = 2)
cosZones <- abs(cospi(Zones/180))
SunWt <- Func(Zones)
Rin <- Incident(S,SunWt)
Sarr <- rep(0,length(Zones))
#Sarr[1] <- Sun7(S,s)
BLACK <- matrix(NA, nrow=length(Zones), ncol=length(Zones))
WHITE <- matrix(NA, nrow=length(Zones), ncol=length(Zones))
for(h in c(1:length(Zones))) {
S <- 920+(h-1)*10
Sarr[h] <- S
T <- gauss(Zones,0,50,31.6)-6
w <- rep(w0,length(Zones)) #0.5
b <- rep(b0,length(Zones)) #0.2
u <- rep(1-w0-b0,length(Zones))
a <- w*aW+b*aB+u*alb(T,ai,ab,gamma,delta)
for(i in c(1:cycles)) {
Rin <- Incident(S,SunWt)
Tcos <- cosZones*T
Tm <- sum(Tcos)/sum(cosZones)
T <- (Rin*(1-a)+K*Tm-A) / (B+K)
Tw <- T+c*(a-aW)
Tb <- T+c*(a-aB)
Fw <- 1-k*(T0-Tw)^2
Fb <- 1-k*(T0-Tb)^2
for(j in c(1:length(Zones))){
if(Fw[j]<0){Fw[j]=0}
if(Fb[j]<0){Fb[j]=0} }
w <- w+w*(u*Fw-D)
b <- b+b*(u*Fb-D)
for(j in c(1:length(Zones))){
if(w[j]<0.001){w[j]=0.001}
if(b[j]<0.001){b[j]=0.001} }
u <- 1-w-b
a <- w*aW+b*aB+u*alb(T,ai,ab,gamma,delta)
}
BLACK[h,] <- b
WHITE[h,] <- w
}
return ( plot_ly( x=Zones, y=Sarr, z=~BLACK ,colors = colorRamp(c("white", "black")), type = "heatmap") )
}
ebm_Db2 <- function(cycles,w0,b0,c,k,D,A,B,K,ai,ab,aW,aB,gamma,delta){
Incident <- function(x,y){ x*y/4 }
Func <- function(x){ 0.7768699*cos(0.0164348*x)^2+0.4617747 }
Sun6 <- function(x){1370*(1+0.1*cospi(x/180))}
Zones <- seq(-89, 89, by = 2)
cosZones <- abs(cospi(Zones/180))
SunWt <- Func(Zones)
Rin <- Incident(S,SunWt)
Sun6 <- function(x){1370*(1+0.1*cospi(x/180))}
T <- gauss(Zones,0,50,31.6)-6
w <- rep(w0,length(Zones)) #0.5
b <- rep(b0,length(Zones)) #0.2
u <- rep(1-w0-b0,length(Zones))
a <- w*aW+b*aB+u*alb(T,ai,ab,gamma,delta)
Barr <- rep(0,cycles)
Warr <- rep(0,cycles)
Uarr <- rep(0,cycles)
Tarr <- rep(0,cycles)
I <- rep(0,cycles)
TEMP <- matrix(NA, nrow=length(Zones), ncol=cycles)
BLACK <- matrix(NA, nrow=length(Zones), ncol=cycles)
WHITE <- matrix(NA, nrow=length(Zones), ncol=cycles)
for(i in c(1:cycles)) {
S <- Sun6(i) # oppure costante S <- 1370
Rin <- Incident(S,SunWt)
Tcos <- cosZones*T
Tm <- sum(Tcos)/sum(cosZones)
T <- (Rin*(1-a)+K*Tm-A) / (B+K)
Tw <- T+c*(a-aW)
Tb <- T+c*(a-aB)
Fw <- 1-k*(T0-Tw)^2
Fb <- 1-k*(T0-Tb)^2
for(j in c(1:length(Zones))){
if(Fw[j]<0){Fw[j]=0}
if(Fb[j]<0){Fb[j]=0} }
w <- w+w*(u*Fw-D)
b <- b+b*(u*Fb-D)
for(j in c(1:length(Zones))){
if(w[j]<0.001){w[j]=0.001}
if(b[j]<0.001){b[j]=0.001} }
u <- 1-w-b
a <- w*aW+b*aB+u*alb(T,ai,ab,gamma,delta)
Barr[i] <- b[45]
Warr[i] <- w[45]
Uarr[i] <- u[45]
I[i] <- i
Tarr[i] <- T[45]
TEMP[,i] <- T
BLACK[,i] <- b
WHITE[,i] <- w
}
return ( BLACK )
}
ebm_Dw1 <- function(cycles,s,w0,b0,c,k,D,A,B,K,ai,ab,aW,aB,gamma,delta){
Incident <- function(x,y){ x*y/4 }
Func <- function(x){ 0.7768699*cos(0.0164348*x)^2+0.4617747 }
Sun7 <- function(x,y){x*y}
Zones <- seq(-89, 89, by = 2)
cosZones <- abs(cospi(Zones/180))
SunWt <- Func(Zones)
Rin <- Incident(S,SunWt)
Sarr <- rep(0,length(Zones))
#Sarr[1] <- Sun7(S,s)
BLACK <- matrix(NA, nrow=length(Zones), ncol=length(Zones))
WHITE <- matrix(NA, nrow=length(Zones), ncol=length(Zones))
for(h in c(1:length(Zones))) {
S <- 920+(h-1)*10
Sarr[h] <- S
T <- gauss(Zones,0,50,31.6)-6
w <- rep(w0,length(Zones)) #0.5
b <- rep(b0,length(Zones)) #0.2
u <- rep(1-w0-b0,length(Zones))
a <- w*aW+b*aB+u*alb(T,ai,ab,gamma,delta)
for(i in c(1:cycles)) {
Rin <- Incident(S,SunWt)
Tcos <- cosZones*T
Tm <- sum(Tcos)/sum(cosZones)
T <- (Rin*(1-a)+K*Tm-A) / (B+K)
Tw <- T+c*(a-aW)
Tb <- T+c*(a-aB)
Fw <- 1-k*(T0-Tw)^2
Fb <- 1-k*(T0-Tb)^2
for(j in c(1:length(Zones))){
if(Fw[j]<0){Fw[j]=0}
if(Fb[j]<0){Fb[j]=0} }
w <- w+w*(u*Fw-D)
b <- b+b*(u*Fb-D)
for(j in c(1:length(Zones))){
if(w[j]<0.001){w[j]=0.001}
if(b[j]<0.001){b[j]=0.001} }
u <- 1-w-b
a <- w*aW+b*aB+u*alb(T,ai,ab,gamma,delta)
}
BLACK[h,] <- b
WHITE[h,] <- w
}
return ( plot_ly( x=Zones, y=Sarr, z=~WHITE ,colors = colorRamp(c("white", "black")), type = "heatmap") )
}
ebm_Dw2 <- function(cycles,w0,b0,c,k,D,A,B,K,ai,ab,aW,aB,gamma,delta){
Incident <- function(x,y){ x*y/4 }
Func <- function(x){ 0.7768699*cos(0.0164348*x)^2+0.4617747 }
Sun6 <- function(x){1370*(1+0.1*cospi(x/180))}
Zones <- seq(-89, 89, by = 2)
cosZones <- abs(cospi(Zones/180))
SunWt <- Func(Zones)
Rin <- Incident(S,SunWt)
Sun6 <- function(x){1370*(1+0.1*cospi(x/180))}
T <- gauss(Zones,0,50,31.6)-6
w <- rep(w0,length(Zones)) #0.5
b <- rep(b0,length(Zones)) #0.2
u <- rep(1-w0-b0,length(Zones))
a <- w*aW+b*aB+u*alb(T,ai,ab,gamma,delta)
Barr <- rep(0,cycles)
Warr <- rep(0,cycles)
Uarr <- rep(0,cycles)
Tarr <- rep(0,cycles)
I <- rep(0,cycles)
TEMP <- matrix(NA, nrow=length(Zones), ncol=cycles)
BLACK <- matrix(NA, nrow=length(Zones), ncol=cycles)
WHITE <- matrix(NA, nrow=length(Zones), ncol=cycles)
for(i in c(1:cycles)) {
S <- Sun6(i) # oppure costante S <- 1370
Rin <- Incident(S,SunWt)
Tcos <- cosZones*T
Tm <- sum(Tcos)/sum(cosZones)
T <- (Rin*(1-a)+K*Tm-A) / (B+K)
Tw <- T+c*(a-aW)
Tb <- T+c*(a-aB)
Fw <- 1-k*(T0-Tw)^2
Fb <- 1-k*(T0-Tb)^2
for(j in c(1:length(Zones))){
if(Fw[j]<0){Fw[j]=0}
if(Fb[j]<0){Fb[j]=0} }
w <- w+w*(u*Fw-D)
b <- b+b*(u*Fb-D)
for(j in c(1:length(Zones))){
if(w[j]<0.001){w[j]=0.001}
if(b[j]<0.001){b[j]=0.001} }
u <- 1-w-b
a <- w*aW+b*aB+u*alb(T,ai,ab,gamma,delta)
Barr[i] <- b[45]
Warr[i] <- w[45]
Uarr[i] <- u[45]
I[i] <- i
Tarr[i] <- T[45]
TEMP[,i] <- T
BLACK[,i] <- b
WHITE[,i] <- w
}
return ( WHITE )
}Solar Luminosity
\[ S_1(t)= \frac{S}{100}t \] \[ S_2(t)= S(\sin^2(t+90)) \] \[ S_3(t)= S\Big(1-\frac{1}{\sqrt{2\pi}}e^{-(t-50)^2/2}\Big) \] \[ S_4(t)=S\Big(1-\frac{1}{3}\delta(t-50)\Big) \] \[ S_5(t)=\frac{1}{100}(|t-150|+25) \] \[ S_6(t)=S\Big(1+\frac{1}{10} cos(t)\Big) \]
Albedo
\[ a(T)=\frac{e^{\gamma(T+\delta)}}{e^{\gamma(T+\delta)}+1} (\beta-\alpha) + \alpha \]
EBM
\[ T =\frac{ R_{in}(1-a(T))+KT_m-A }{ B+K } \]
Daisyworld
\[ a(T)=wa_w+ba_b+u\Big(\frac{e^{\gamma(T+\delta)}}{e^{\gamma(T+\delta)}+1} (\beta-\alpha) + \alpha\Big)\] \[ T_w=T+c(a-a_w) \] \[ T_b=T+c(a-a_b) \] \[ F_w=1-k(T_0-T_w)^2 \] \[ F_b=1-k(T_0-T_b)^2 \] \[ w'=w+w(uF_w-D) \] \[ b'=b+b(uF_b-D) \]
plot_albedo <- ggplot(data.frame(x= seq(-15,15, by=0.1),y=alb(seq(-15,15, by=0.1),ai,ab,gamma,delta)))+geom_line(aes(x,y),colour="blue")+xlab("Temperature")+ylab("Albedo")
plot_albedo
Run 0
data01 <- ebm01(150,A,B,K,ai,ab,gamma,delta) #from 300
data02 <-ebm02(150,A,B,K,ai,ab,gamma,delta) #from 300
plot01 <- ggplot(data01,aes(Zones)) +
geom_line(aes(y=T, colour = "Temperature",linetype="Temperature")) +
geom_line(aes(y=a*25, colour = "Albedo",linetype="Albedo"))+
geom_line(aes(y=Ti, colour = "Initial Temperature",linetype="Initial Temperature"))+xlab("Latitude")+ylab("Temperature")+scale_colour_manual(name="Legend",values=c("blue","red","red"))+scale_linetype_manual(name="Legend",values=c("Temperature"=1, "Albedo"=1, "Initial Temperature"=2))
plot02 <- ggplot(data02) +
geom_line(aes(t, Temperature),colour = 'green')+xlab("Time")+ylab("Mean Temperature")
plot01
plot02
Run 1
data11 <- ebm11(100,300,A,B,K,ai,ab,gamma,delta)
plot11 <- ggplot(data11,aes(J,Temp))+geom_line(aes(J, Temp),colour = 'red')+xlab("Time")+ylab("Mean Temperature")+ggtitle("Non Dynamic Mean Temperature (T(t) independent from T(t-1))")
plot11
TEMP1 <- ebm12(100,300,A,B,K,ai,ab,gamma,delta)
plot_ly(z=~TEMP1)%>% add_surface() %>% layout(
title = "With Initialization", scene = list(
xaxis = list(title = "S/100"),
yaxis = list(title = "Latitude"),
zaxis = list(title = "T") )) Run 2
data21 <- ebm21(180,300,A,B,K,ai,ab,gamma,delta )
plot21 <- ggplot(data21,aes(J,Temp2))+geom_line(aes(J, Temp2),colour = 'red')+xlab("Time")+ylab("Mean Temperature")+ggtitle("Sinusoidal Perturbation")
plot21
TEMP2 <- ebm22(180,300,A,B,K,ai,ab,gamma,delta )
plot_ly(z=~TEMP2) %>% add_surface() %>% layout(
title = "Without Initialization",scene = list(
xaxis = list(title = "S/100"),
yaxis = list(title = "Latitude"),
zaxis = list(title = "T") ))Run 3 & 4
data31 <- ebm31(200,300,A,B,K,ai,ab,gamma,delta )
data41 <- ebm41(200,300,A,B,K,ai,ab,gamma,delta )
plot31 <- ggplot(data31,aes(J,Temp))+geom_line(aes(J, Temp),colour = 'red')+xlab("Time")+ylab("Mean Temperature")+ggtitle("Gaussian Perturbation")
plot41 <- ggplot(data41,aes(J,Temp))+geom_line(aes(J, Temp),colour = 'red')+xlab("Time")+ylab("Mean Temperature")+ggtitle("Delta Perturbation")
plot31
plot41
TEMP3 <- ebm32(200,300,A,B,K,ai,ab,gamma,delta )
TEMP4 <- ebm42(200,300,A,B,K,ai,ab,gamma,delta)
plot_ly(z=~TEMP3) %>% add_surface()%>% layout(
title = "Without Initialization", scene = list(
xaxis = list(title = "S/100"),
yaxis = list(title = "Latitude"),
zaxis = list(title = "T") ))plot_ly(z=~TEMP4) %>% add_surface()%>% layout(
title = "Without Initialization",scene = list(
xaxis = list(title = "S/100"),
yaxis = list(title = "Latitude"),
zaxis = list(title = "T") ))Hysteresis Cycles
data23 <- ebm23(180,300,A,B,K,ai,ab,gamma,delta )
plot23 <- ggplot(data23,aes(Sarr,Temp2))+geom_point(aes(Sarr, Temp2),colour = 'yellow')+ylab("Mean Temperature")+xlab("S_2(t)")+ggtitle("Mean temperature vs. Sinusoidally Fluctuating Incoming Radiation")
plot23
data51 <- ebm51(300,60,A,B,K,ai,ab,gamma,delta )
plot51 <- ggplot(data51,aes(Sarr,Temp5,group=1))+geom_point(aes(Sarr, Temp5),colour = 'yellow')+geom_segment(aes(x = Sarr[89], y = Temp5[89], xend = Sarr[90], yend = Temp5[90]),colour = 'yellow')+geom_segment(aes(x = Sarr[258], y = Temp5[258], xend = Sarr[259], yend = Temp5[259]),colour = 'yellow')+ylab("Mean Temperature")+xlab("S_5(t)")+ggtitle("Mean temperature vs. Linearly Fluctuating Incoming Radiation")
plot51
Embedding Daisies in EBM
# WITHOUT DAISIES
data_ND1 <- ebm_ND1(150,w0,b0,A,B,K,ai,ab,aW,aB,gamma,delta) #from 500
data_ND2 <- ebm_ND2(150,w0,b0,A,B,K,ai,ab,aW,aB,gamma,delta)
TEMP <- ebm_ND3(150,w0,b0,A,B,K,ai,ab,aW,aB,gamma,delta)
plot_ND1 <- ggplot(data_ND1,aes(Zones))+geom_line(aes(y=w, colour = "% White"))+geom_line(aes(y=b, colour = "% Black"))+geom_line(aes(y=u, colour="% Bare Ground"))+geom_line(aes(y=T/5, colour="Temperature"))+ylab("Temperature & Albedo")+xlab("Latitude")+ggtitle("Without Daisies")+scale_colour_manual(name="Legend",values=c("brown","black","white", "red"))
plot_ND2 <- ggplot(data_ND2,aes(I))+geom_line(aes(y=Barr,colour="% Black"))+ geom_line(aes(y=Warr,colour="% White"))+ geom_line(aes(y=Uarr, colour="% Bare Ground"))+geom_line(aes(y=Tarr/25,colour="Temperature"))+ylab("Temperature & Albedo")+xlab("Time")+ggtitle("Without Daisies")+scale_colour_manual(name="Legend",values=c("brown","black","white", "red"))
plot_ND1 
plot_ND2
plot_ly(z=~TEMP)%>% add_surface() %>% layout( title="Temperature Without Daisies",
scene = list(
xaxis = list(title = "Time"),
yaxis = list(title = "Latitude"),
zaxis = list(title = "T") )) #WITH DAISIES
data_D1 <- ebm_D1(150,w0,b0,c,k,D,A,B,K,ai,ab,aW,aB,gamma,delta)
data_D2 <- ebm_D2(150,w0,b0,c,k,D,A,B,K,ai,ab,aW,aB,gamma,delta)
TEMP <- ebm_D3(150,w0,b0,c,k,D,A,B,K,ai,ab,aW,aB,gamma,delta)
BLACK <- ebm_Db2(150,w0,b0,c,k,D,A,B,K,ai,ab,aW,aB,gamma,delta)
WHITE <- ebm_Dw2(150,w0,b0,c,k,D,A,B,K,ai,ab,aW,aB,gamma,delta)
plot_D1 <- ggplot(data_D1 ,aes(Zones))+geom_line(aes(y=w, colour = "% White"))+geom_line(aes(y=b, colour = "% Black"))+geom_line(aes(y=u, colour="% Bare Ground"))+geom_line(aes(y=T/5, colour="Temperature"))+ylab("Temperature & Albedo")+xlab("Latitude")+ggtitle("With Daisies")+scale_colour_manual(name="Legend",values=c("brown","black","white", "red"))
plot_D2 <- ggplot(data_D2,aes(I))+geom_line(aes(y=Barr,colour="% Black"))+ geom_line(aes(y=Warr,colour="% White"))+ geom_line(aes(y=Uarr, colour="% Bare Ground"))+geom_line(aes(y=Tarr/25,colour="Temperature"))+ylab("Temperature & Albedo")+xlab("Time")+ggtitle("With Daisies")+scale_colour_manual(name="Legend",values=c("brown","black","white", "red"))
plot_D1
plot_D2
plot_ly(z=~TEMP)%>% add_surface() %>% layout(title="Temperature With Daisies",
scene = list(
xaxis = list(title = "Time"),
yaxis = list(title = "Latitude"),
zaxis = list(title = "T") ))plot_ly(z=~BLACK)%>% add_surface() %>% layout(title="Distribution of Black Daisies",
scene = list(
xaxis = list(title = "Time"),
yaxis = list(title = "Latitude"),
zaxis = list(title = "Black") )) plot_ly(z=~WHITE)%>% add_surface() %>% layout(title="Distribution of White Daisies",
scene = list(
xaxis = list(title = "Time"),
yaxis = list(title = "Latitude"),
zaxis = list(title = "White") ))ebm_Db1(150,1370/920,w0,b0,c,k,D,A,B,K,ai,ab,aW,aB,gamma,delta)%>% layout(title="Distribution of Black Daisies",
xaxis = list(title = "Latitude"),
yaxis = list(title = "Solar Luminosity") )ebm_Dw1(150,1370/920,w0,b0,c,k,D,A,B,K,ai,ab,aW,aB,gamma,delta) %>% layout(title="Distribution of White Daisies",
xaxis = list(title = "Latitude"),
yaxis = list(title = "Solar Luminosity") )